##########################################################################################

library(dplyr)
library(ggplot2)
library(data.table)
library(RColorBrewer)
library(optparse)

##########################################################################################

option_list <- list(
    make_option(c("--maf_cancer_file"), type = "character") ,
    make_option(c("--maf_im_file"), type = "character") ,
    make_option(c("--images_path"), type = "character") ,
    make_option(c("--info_file"), type = "character")
)

if(1!=1){
    maf_cancer_file <- "~/20220915_gastric_multiple/dna_combinePublic/maf_public/All_use.maf"
    maf_im_file <- "~/20220915_gastric_multiple/dna_combinePublic/maf_public/All_use.IM.maf"
    info_file <- "~/20220915_gastric_multiple/dna_combinePublic/public_ref/combine/MutationInfo.combine.addMolecularSubType.rmMIX.tsv"
}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

maf_cancer_file <- opt$maf_cancer_file
maf_im_file <- opt$maf_im_file
info_file <- opt$info_file
images_path <- opt$images_path

dir.create(images_path , recursive = T)

###########################################################################################

info <- data.frame(fread(info_file))
dat_maf_cancer <- data.frame(fread( maf_cancer_file ))
dat_maf_im <- data.frame(fread( maf_im_file ))

###########################################################################################

Variant_Type <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")

###########################################################################################

info <- subset( info , !Age %in% c("[Not Available]" , "unknown") )
info$Age <- as.numeric(info$Age)
info$Age_divide <- ifelse(info$Age > 66 , "Older" , "Younger")

###########################################################################################
## 不看IGC + DGC 均有的样本
## 只看MSS的样本
info_gc_divide <- subset( info , Class != "IGC + DGC" & MS_Type == "MSS" ) 
info_gc <- info_gc_divide
info_gc$Class <- "GC"
info_im <- subset( info , Class != "IGC + DGC" & From == "NJMU" & MS_Type == "MSS" ) 
info_im$Class <- "IM"

info_use <- rbind( info_gc_divide , info_gc , info_im )
info_all <- info_use
info_all$From <- 'All'
info_use <- rbind(info_all , info_use)

###########################################################################################
## 注释肿瘤样本的类型
maf <- dat_maf_cancer
maf_use <- data.frame(
    Hugo_Symbol = maf$Hugo_Symbol, 
    Chromosome = maf$Chromosome , Start_Position =  maf$Start_position , 
    Reference_Allele = maf$Reference_Allele , Tumor_Seq_Allele2 = maf$Tumor_Seq_Allele2 , 
    Variant_Classification = maf$Variant_Classification , Tumor = maf$Tumor_Sample_Barcode , From = maf$From )

maf_use <- subset( maf_use , Variant_Classification %in% Variant_Type )

## 标记年龄分组
maf_use <- merge( maf_use , info[,c("Tumor" , "Class" , "Age_divide" )] , by = "Tumor")
maf_cancer <- maf_use

###########################################################################################
## 肠化样本
maf <- dat_maf_im
maf_use <- data.frame(
    Hugo_Symbol = maf$Hugo_Symbol, 
    Chromosome = maf$Chromosome , Start_Position =  maf$Start_position , 
    Reference_Allele = maf$Reference_Allele , Tumor_Seq_Allele2 = maf$Tumor_Seq_Allele2 , 
    Variant_Classification = maf$Variant_Classification , Tumor = maf$Tumor_Sample_Barcode , From = maf$From )

maf_use <- subset( maf_use , Variant_Classification %in% Variant_Type )
maf_im <- maf_use
maf_im$Class <- "IM"
maf_im <- subset( maf_im , Tumor %in% info_im$Tumor  )
maf_im <- merge( maf_im , info[,c("Tumor" , "Age_divide")] , by = "Tumor")

###########################################################################################
## 总的，不同来源的
## IGC和DGC分开的，合并的
maf_cancer2 <- maf_cancer
maf_cancer2$Class <- "GC"

maf_use <- rbind(maf_cancer , maf_cancer2 , maf_im)
maf_use$From <- "All"
maf_use <- rbind(maf_use , maf_cancer , maf_cancer2 , maf_im)

maf_use <- subset( maf_use , Tumor %in% info_use$Tumor  )

###########################################################################################
## 按照病理类型计算突变率
## 以人为单位
mut_rate <- maf_use %>%
group_by( Hugo_Symbol , Class , Age_divide ) %>%
summarize( MutNum = length(unique(Tumor )) )
mut_rate$id <- paste0( mut_rate$Class , ":" , mut_rate$Age_divide )

class_num <- info_use %>%
group_by( Class , Age_divide) %>%
summarize( SampleNum = length(unique(Tumor)) )
class_num$id <- paste0( class_num$Class , ":" , class_num$Age_divide )

mut_rate <- merge( mut_rate , class_num[,c("id" , "SampleNum")] , by = "id" )
## 分IGC、DGC和IM
mut_rate$MutRate <- mut_rate$MutNum/mut_rate$SampleNum

out_name <- paste0( images_path , "/MutRate.Age_divide.tsv" )
write.table( mut_rate , out_name , row.names = F , quote = F , sep = "\t" )

###########################################################################################
## 计算关注的驱动基因的突变数量
gene_order <- c("TP53","ARID1A","CDH1","APC","SMAD4","MUC6","PIK3CA",
              "CTNNB1","RHOA","ERBB2","CFTR","KRAS","MAP2K7","ARID2",
              "RNF43","TGFBR2","BMP6","FBXW7","CDKN2A","MTRR")

maf_dirver_use <- subset(maf_use , Hugo_Symbol %in% gene_order)

dat_info_combine2 <- maf_dirver_use %>%
group_by( Tumor , Class , Age_divide ) %>%
summarize( MutNum=length(unique(Start_Position)) )

## 画图
trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("P == ",down," %*% 10","^",up)
    return(text)
}

dat_tmp <- c()
for( class in unique(dat_info_combine2$Class) ){

    dat <- data.frame(subset( dat_info_combine2 , Class == class ))

    a <- dat[dat$Age_divide==unique(dat$Age_divide)[1],"MutNum"]
    b <- dat[dat$Age_divide==unique(dat$Age_divide)[2],"MutNum"]
    
    p <- wilcox.test( a , b )$p.value
    if( p < 0.01 ){
        p_text <- trans(p)
    }else{
        p_text <- paste0( "P == " , round(as.numeric(p) , 3) ) 
    }
    dat$p_text <- ""
    dat$p_text[1] <- p_text
    dat_tmp <- rbind(dat_tmp , dat)
}

dat_tmp <- subset( dat_tmp , Class != "GC" )
dat_tmp$Class <- factor( dat_tmp$Class , levels = c("IM" , "IGC" , "DGC") , order = T )

###########################################################################################

y_lab <- "Number of mutations in SMGs"
y_max <- max(dat_tmp$MutNum)
col <- c(
    rgb(red=247,green=184,blue=71,alpha=255,max=255) ,
    rgb(red=2,green=100,blue=190,alpha=255,max=255) 
    )

###########################################################################################

plot <- ggplot(dat_tmp, aes(x = Age_divide , y = MutNum, color = Age_divide)) +
    #geom_violin(trim = FALSE) +
    geom_boxplot(size = 1.2 , outlier.alpha=0) + ## 去除散点，加粗线
    geom_jitter(position = position_jitterdodge(0.8) , size = 1) + 
    scale_color_manual(values=col) +
    facet_grid( .~ Class , scales = "free_x" ) +
    xlab(NULL) +
    ylab(y_lab)+
    theme_bw() +
    geom_text(aes(label=p_text , y = y_max , x = 1.5 ) , parse = TRUE,size=2 , color = "black") +
    #geom_text(aes(label=p_text , y = y_max),size=5 , color = "black") +
    theme(
        legend.position = 'none',
        legend.title = element_blank() ,
        panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        panel.background = element_blank(),
        panel.border = element_blank(),
        plot.title = element_text(size = 6,color="black",face='bold'),
        legend.text = element_text(size = 3,color="black",face='bold'),
        axis.text.y = element_text(size = 6,color="black",face='bold'),
        axis.title.x = element_text(size = 6,color="black",face='bold'),
        axis.title.y = element_text(size = 6,color="black",face='bold'),
        axis.text.x = element_text(size = 6,color="black",face='bold'),
        axis.line = element_line(size = 0.5)) 

out_name <- paste0(images_path , "/DriverMutCompare.Age_divide.pdf")
ggsave(file=out_name,plot=plot,width=1.79/1.3,height=2.34/1.3)

########################################################################
## 柱状图
## 中位数和IQR
ebtop<-function(x){
  return(quantile(x)[4])
}
ebbottom<-function(x){
  return(quantile(x)[2])
}

y_max <- max(dat_tmp$MutNum)

plot <- ggplot(dat_tmp, aes(x = Age_divide , y = MutNum, fill = Age_divide)) +
scale_fill_manual(values=col) +
scale_color_manual(values=col) +
stat_summary(geom = "bar",fun = "median", color = "black" ,
           position = position_dodge(0.9))+
stat_summary(geom = "errorbar",
           fun.min = ebbottom,
           fun.max = ebtop,
           position = position_dodge(0.9),
           width=0.2)+
scale_y_continuous(expand = expansion(mult = c(0,0.1)))+
facet_grid( .~ Class , scales = "free_x" ) +
geom_text(aes(label=p_text , y = y_max + 0.05 , x = 1.5) , size=5 , color = "black" , parse = TRUE) +
theme_bw()+
xlab(NULL) +
ylab(y_lab)+
ylim(-0.02 , y_max + 0.07) +
theme(panel.grid = element_blank())+
theme(
    legend.position = "none",
    legend.title = element_blank() ,
    panel.grid.major=element_blank(),
    panel.grid.minor=element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(),
    plot.title = element_text(size = 12,color="black",face='bold'),
    legend.text = element_text(size = 8,color="black",face='bold'),
    axis.text.y = element_text(size = 8,color="black",face='bold'),
    axis.title.x = element_text(size = 10,color="black",face='bold'),
    axis.title.y = element_text(size = 10,color="black",face='bold'),
    axis.text.x = element_text(size = 10,color="black",face='bold') ,
    axis.line = element_line(size = 0.5)) 

out_name <- paste0(images_path , "/DriverMutCompare.Age_divide.pdf")
ggsave(file=out_name,plot=plot,width=4,height=3)

dat_tmp %>%
group_by(Age_divide , Class) %>%
summarize(d25 = quantile(MutNum)[1] , d50 = quantile(MutNum)[2] , d75 = quantile(MutNum)[3])
